2.2.2 Data exploration

summaries <- learning2014_csv %>%
  summarise(across(where(is.numeric), list(mean = ~mean(.), sd = ~sd(.), min = ~min(.), max = ~max(.), median = ~median(.))))

print(t(summaries))

##                       [,1]
## age_mean        25.5120482
## age_sd           7.7660785
## age_min         17.0000000
## age_max         55.0000000
## age_median      22.0000000
## attitude_mean    3.1427711
## attitude_sd      0.7299079
## attitude_min     1.4000000
## attitude_max     5.0000000
## attitude_median  3.2000000
## deep_mean        3.6797189
## deep_sd          0.5541369
## deep_min         1.5833333
## deep_max         4.9166667
## deep_median      3.6666667
## stra_mean        3.1212349
## stra_sd          0.7718318
## stra_min         1.2500000
## stra_max         5.0000000
## stra_median      3.1875000
## surf_mean        2.7871486
## surf_sd          0.5288405
## surf_min         1.5833333
## surf_max         4.3333333
## surf_median      2.8333333
## points_mean     22.7168675
## points_sd        5.8948836
## points_min       7.0000000
## points_max      33.0000000
## points_median   23.0000000

ggpairs(learning2014_csv, mapping = aes(), lower = list(combo = wrap("facethist", bins = 20))) + theme_minimal()

2.2.3 Linear regression

I have chosen the three variables to check: attitude, age and surf

linear_modelling_for_learning2014_csv  <- lm(points ~ surf + attitude + age, data = learning2014_csv) 
summary(linear_modelling_for_learning2014_csv)

## 
## Call:
## lm(formula = points ~ surf + attitude + age, data = learning2014_csv)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.973  -3.430   0.167   3.997  10.444 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 16.85628    3.53786   4.765 4.18e-06 ***
## surf        -0.96049    0.79943  -1.201    0.231    
## attitude     3.42388    0.57353   5.970 1.45e-08 ***
## age         -0.08713    0.05361  -1.625    0.106    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.294 on 162 degrees of freedom
## Multiple R-squared:  0.2082, Adjusted R-squared:  0.1935 
## F-statistic:  14.2 on 3 and 162 DF,  p-value: 2.921e-08

Only one explanatory variable showed statistically significant connection to the points outcome, thus only attitude variable will be kept for the next modelling step.

linear_modelling_for_learning2014_csv_pruned  <-lm(points ~ attitude, data = learning2014_csv) 
summary(linear_modelling_for_learning2014_csv_pruned)

## 
## Call:
## lm(formula = points ~ attitude, data = learning2014_csv)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.9763  -3.2119   0.4339   4.1534  10.6645 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.6372     1.8303   6.358 1.95e-09 ***
## attitude      3.5255     0.5674   6.214 4.12e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.32 on 164 degrees of freedom
## Multiple R-squared:  0.1906, Adjusted R-squared:  0.1856 
## F-statistic: 38.61 on 1 and 164 DF,  p-value: 4.119e-09

The second model is simpler, but have not lost too much of the explanatory power: - The p-value is comparably low, much lower than 0.001, meaning that there is strong support to reject the hypothesis that the 2 variables have no connection. - The estimate shows how the two variables are connected, for one unit change in points the attitude variable changes by 3.5255 The adjusted R-value 0.1856 is relatively low, indicating that the model may not be very effective in predicting or explaining the dependent variable. This might be due to various reasons like missing important predictors, non-linear relationships, or high levels of noise in the data.

2.2.3.1 Linear model diagnositic plots

plot(linear_modelling_for_learning2014_csv_pruned, which=c(1,2,5))

Using these plots we can investigate the assumptions of the model.

1. We can use the “Residuals vs fitted” plots to investigate the constant variance assumption. If the variance is constant, we should expect to see point distributed without noticeable structure i.e. randomly. This is more or less what we see in our data. Although the points around 28 on the x axis seem to be bunched up, this may also be the result of low n-numbers, so we can assume constant variance.
2. The q-q recapitulates plot can be used to identify whether the “normal distribution of residuals” is met. The middle part of the plot follows the 45 degree line, meaning that the distribution is close to normal for the bulk of the data. The tails of the plot do deviate form the line, meaning that there is a deviation form the perfect normality. At the bottom left there are 3 points that can be considered outliers. Further investigation is needed to determine whether these point should not be included into furhter analysis.
3. The residuals vs Leverage plot shows whitewater the data contains influential outliers. Three points are highlighted: 71, 56 and 35. These should be investigated to determine if the data is fine, for example there can be a data entry mistake, missing value problem, or a true outlier. The points 56 and 35 are present on both the q-q and the residuals vs leverage plots, so the next step would be checking the data to see what happened to these points.

3.3 Choosing the variable

I decided to look at “failures”, “absences”, “freetime” and “famrel”. I would expect that:

1. The failures might be positively associated with the higher consumption
2. The abscences might also be positively associated with the higher consumption
3. The freetime may be associated with weekly consuption, but It is possible that the association will be inverse.
4. The family relations may be inversely correlated.

3.4 Testing the assumptions:

summaries <- aa_alc %>% group_by(high_use) %>% select(.,c("failures", "absences", "freetime", "famrel","high_use")) %>% 
  summarise(across(where(is.numeric), list(Average = ~mean(.), sd = ~sd(.))))

gt(summaries) %>% cols_align("left") %>% opt_stylize(color="gray", style=3)  %>% tab_header("Table for chosen variables")

These are the averages and standard deviation values for all the variables, to better see which variables are interesting to look at closer.

First looking at overall structure, whehther there is something to connect the chosen variables.

ggpairs(select(aa_alc,c("failures", "absences", "freetime", "famrel","high_use")), mapping = aes(), lower = list(combo = wrap("facethist", bins = 20))) + theme_minimal()

As it turns out the free time is strongly correlated with family relations, which I did not expect.

ggpairs(select(aa_alc,c("failures","high_use")), mapping = aes(), lower = list(combo = wrap("facethist", bins = 20))) + theme_minimal()

ggpairs(select(aa_alc,c("absences","high_use")), mapping = aes(), lower = list(combo = wrap("facethist", bins = 20))) + theme_minimal()

ggpairs(select(aa_alc,c("freetime","high_use")), mapping = aes(), lower = list(combo = wrap("facethist", bins = 20))) + theme_minimal()

ggpairs(select(aa_alc,c("famrel","high_use")), mapping = aes(), lower = list(combo = wrap("facethist", bins = 20))) + theme_minimal()

3.5 Using logistic regression

first_log_regression <- glm(high_use ~ failures + absences + freetime + famrel, data = aa_alc, family = "binomial")

# create model summary and extract the coeffs
summary(first_log_regression)

## 
## Call:
## glm(formula = high_use ~ failures + absences + freetime + famrel, 
##     family = "binomial", data = aa_alc)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.55560    0.63703  -2.442 0.014608 *  
## failures     0.58297    0.20259   2.878 0.004007 ** 
## absences     0.08375    0.02230   3.755 0.000173 ***
## freetime     0.43091    0.12804   3.365 0.000764 ***
## famrel      -0.31981    0.13098  -2.442 0.014622 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 452.04  on 369  degrees of freedom
## Residual deviance: 408.17  on 365  degrees of freedom
## AIC: 418.17
## 
## Number of Fisher Scoring iterations: 4

From the model we can see the following:

• The higher values of failures absences and the freetime predictors are associated with a higher likelihood of high_use being TRUE, for family relations the opposite is true.
• All predictors appear to be statistically significant, as indicated by their p-values and significance codes
• The reduction in deviance from null to residual suggests the model with predictors fits better than the null model.
• The reduction in deviance from null to residual suggests the model with predictors fits better than the null model - the one without ony predictor variables.

coeffs <- summary(first_log_regression) %>% coefficients()
coeffs

##                Estimate Std. Error   z value     Pr(>|z|)
## (Intercept) -1.55560199 0.63703172 -2.441954 0.0146080178
## failures     0.58296988 0.20258843  2.877607 0.0040070412
## absences     0.08375499 0.02230317  3.755295 0.0001731376
## freetime     0.43090906 0.12804024  3.365419 0.0007642749
## famrel      -0.31981102 0.13098388 -2.441606 0.0146221015

OddsRatio <- exp(coeffs)
confidence <- confint(first_log_regression)

## Waiting for profiling to be done...

result <- cbind(OddsRatio, confidence)
result

##              Estimate Std. Error     z value Pr(>|z|)      2.5 %      97.5 %
## (Intercept) 0.2110623   1.890860  0.08699073 1.014715 -2.8317664 -0.32620631
## failures    1.7913506   1.224568 17.77169287 1.004015  0.1917587  0.99258448
## absences    1.0873624   1.022554 42.74681512 1.000173  0.0420769  0.13014453
## freetime    1.5386556   1.136599 28.94562433 1.000765  0.1840339  0.68715409
## famrel      0.7262863   1.139949  0.08702100 1.014730 -0.5789180 -0.06341917

Here we can see the same, the first three predictors have positive association with consumption, the family relation variable has a negative association. We can see that the strongest connection is with failures, then freetime and famrel, while the link to absences is low.

It is important to keep in mind that this is logistical regression, not the previosuly investigated linear regression. Tha means that the estimate in this table represents the “odds ratos” and be thought of in terms of likelihood. It means, that for examples in our case the increase in failures by one unit increases the likelihood of high alcohol consumption by 1.8.

The results are in agreement with the hypotheses I started with. Interestingly there seems to be a connection between the free time and the consumption levels, which were not obvious when just looking at the data.

3.6 Exploring the predictive power

We can use test the predictive of the model that we created earlier to see if it can be used to actually predict the alcohol consumption based on the the four chosen variables. We can predict for each student the probability of increased consumption.

aa_alc$predicted_probabilities <- predict(first_log_regression, type = "response")  

aa_alc <- mutate(aa_alc,prediction  = predicted_probabilities > 0.5) 

table(Actual = aa_alc$high_use, Predicted = aa_alc$prediction)

##        Predicted
## Actual  FALSE TRUE
##   FALSE   236   23
##   TRUE     85   26

These are the 2x2 cross tabulation count results.

((table(Actual_percentage = aa_alc$high_use, Predicted_percentage = aa_alc$prediction))/length(aa_alc$prediction)) * 100

##                  Predicted_percentage
## Actual_percentage     FALSE      TRUE
##             FALSE 63.783784  6.216216
##             TRUE  22.972973  7.027027

The False outcome was correctly predicted for 236 participants (63%), True was predicted for 26 (7%).

library(ggplot2)
ggplot(aa_alc, aes(x=high_use, y=prediction)) +
  geom_jitter(color="blue") +
  theme_minimal() +
  labs(title="Actual vs Predicted Alcohol Consumption")

confusion_matrix <- table(Predicted = aa_alc$prediction,Actual = aa_alc$high_use)
training_error <- 1 - sum(diag(confusion_matrix)) / sum(confusion_matrix)
training_error

## [1] 0.2918919

The total training error was 29%

sum(aa_alc$high_use)/length(aa_alc$high_use)

## [1] 0.3

simple_guess_accuracy <- mean(aa_alc$high_use == F)
model_accuracy <- mean(aa_alc$high_use == aa_alc$prediction)

simple_guess_accuracy

## [1] 0.7

model_accuracy

## [1] 0.7081081

loss_func <- function(class, prob) {
  n_wrong <- abs(class - prob) > 0.5
  mean(n_wrong)
}
library(boot)

cv <- cv.glm(data = aa_alc, cost = loss_func, glmfit = first_log_regression, K = 10)

cv$delta[1]

## [1] 0.2972973

if (!require(caret) == T) {
  install.packages("caret")
}

## Loading required package: caret

## Loading required package: lattice

## 
## Attaching package: 'lattice'

## The following object is masked from 'package:boot':
## 
##     melanoma

## 
## Attaching package: 'caret'

## The following object is masked from 'package:purrr':
## 
##     lift

library(caret)

if (!require(glmnet) == T) {
  install.packages("glmnet")
}

## Loading required package: glmnet

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## Loaded glmnet 4.1-8

library(glmnet)

predictors <- setdiff(colnames(aa_alc), c("Dalc","Walc","average_alc","high_use","predicted_probabilities","prediction" ))
# Initialize vectors to store errors
training_errors <- c()
testing_errors <- c()
number_of_predictors <- c()
model_summaries <- list()

aa_alc$high_use <- as.factor(aa_alc$high_use)

# Loop over models with decreasing number of predictors
for (i in length(predictors):1) {
    # Select a subset of predictors
    current_predictors <- predictors[1:i]
    formula <- as.formula(paste("high_use ~", paste(current_predictors, collapse = "+")))
    set.seed(123)  # for reproducibility
    fitControl <- trainControl(method = "cv", number = 10)  # 10-fold cross-validation
    model <- train(formula, data = aa_alc, method = "glm", family = "binomial", trControl = fitControl)
    
    # Store training and testing errors
    training_errors <- c(training_errors, mean(model$results$Accuracy))
    testing_errors <- c(testing_errors, 1 - max(model$results$Accuracy))# Replace 'Accuracy' with the appropriate metric
    number_of_predictors <- c(number_of_predictors, i)
    model_summaries[[i]] <- model$coefnames
}

results_df <- data.frame(
    Number_of_Predictors = number_of_predictors,
    Training_Error = training_errors,
    Testing_Error = testing_errors
)

# Plotting
ggplot(results_df, aes(x = Number_of_Predictors)) +
    geom_line(aes(y = Training_Error, colour = "Training Error")) +
    geom_line(aes(y = Testing_Error, colour = "Testing Error")) +
    labs(title = "Training and Testing Errors by Number of Predictors",
         x = "Number of Predictors",
         y = "Error") +
    theme_minimal()

setdiff(model_summaries[[26]],model_summaries[[25]])

## [1] "failures"

predictors  <- rev(predictors)
# Initialize vectors to store errors
training_errors <- c()
testing_errors <- c()
number_of_predictors <- c()
model_summaries <- list()

# Loop over models with decreasing number of predictors
for (i in length(predictors):1) {
    # Select a subset of predictors
    current_predictors <- predictors[1:i]
    formula <- as.formula(paste("high_use ~", paste(current_predictors, collapse = "+")))
    set.seed(123)  # for reproducibility
    fitControl <- trainControl(method = "cv", number = 10)  # 10-fold cross-validation
    model <- train(formula, data = aa_alc, method = "glm", family = "binomial", trControl = fitControl)
    
    # Store training and testing errors
    training_errors <- c(training_errors, mean(model$results$Accuracy))
    testing_errors <- c(testing_errors, 1 - max(model$results$Accuracy))# Replace 'Accuracy' with the appropriate metric
    number_of_predictors <- c(number_of_predictors, i)
    model_summaries[[i]] <- model$coefnames
}

results_df <- data.frame(
    Number_of_Predictors = number_of_predictors,
    Training_Error = training_errors,
    Testing_Error = testing_errors
)

# Plotting
ggplot(results_df, aes(x = Number_of_Predictors)) +
    geom_line(aes(y = Training_Error, colour = "Training Error")) +
    geom_line(aes(y = Testing_Error, colour = "Testing Error")) +
    labs(title = "Training and Testing Errors by Number of Predictors",
         x = "Number of Predictors",
         y = "Error") +
    theme_minimal()

## [1] 0.2454599

predictors[20]

## [1] "guardian"

4.1.4.1: Creating a categorical variable from scaled crime rate

Using quantiles as breakpoints, and removing the old crime rate variable.

quantiles <- quantile(scaled_Boston[, "crim"], probs = c(0, 0.25, 0.5, 0.75, 1))
scaled_Boston$categorical_crim <- cut(scaled_Boston[, "crim"], breaks = quantiles, 
                               include.lowest = TRUE, labels = c("low", "med_low", "med_high", "high"))
scaled_Boston <- scaled_Boston[,-which(names(Boston) == "crim")]

Splitting the Dataset into Train and Test Sets

80% of the data is now in the training set, and the remaining 20% is in the test set.

set.seed(123) # For reproducibility
indexes <- createDataPartition(scaled_Boston$categorical_crim, p = 0.8, list = FALSE)
train_set <- scaled_Boston[indexes, ]
test_set <- scaled_Boston[-indexes, ]

5.2.1 Graphical overview and description of the data

human_data_ <- read.csv(file="/home/alex/ods_course/IODS-project_R/data/human.csv")
human_data_ <- column_to_rownames(human_data_,var ="Country") 

Graphical summary:

ggpairs(human_data_, progress = F)

1. All the distributions are skewed, meaning that extreme values are more common for some variables. Life expectancy tends to high, GNI_C tend to low. 2. A lot of variables are strongly correlated positively and negatively: - Life expectancy, education levels, GNI are positively correlated with good significance, and are negatively correlated with “Maternal Mortality Ratio” and “Adolescent Birth Rate” with similar significance. This is to be expected. - MMortality is correlated with adolscent births, which is to be expected. - Higher secondary education in females is positivly correlated with GNI and life expectancy.

These resutls are not surprising, as a lot of the indexes are associated with the overall wealth, even though there might be outliers.

Summary:

summary(human_data_)

##   FM_2ED_Ratio    FM_LFPR_Ratio        EduExp        LifeExp_B    
##  Min.   :0.1717   Min.   :0.1857   Min.   : 5.40   Min.   :49.00  
##  1st Qu.:0.7254   1st Qu.:0.5981   1st Qu.:11.22   1st Qu.:66.42  
##  Median :0.9327   Median :0.7514   Median :13.45   Median :74.00  
##  Mean   :0.8499   Mean   :0.7034   Mean   :13.13   Mean   :71.58  
##  3rd Qu.:0.9958   3rd Qu.:0.8525   3rd Qu.:15.20   3rd Qu.:76.95  
##  Max.   :1.4967   Max.   :1.0380   Max.   :20.20   Max.   :83.50  
##      GNI_C           MMRatio          AdBRate            RP_p      
##  Min.   :   581   Min.   :   1.0   Min.   :  0.60   Min.   : 0.00  
##  1st Qu.:  4495   1st Qu.:  13.0   1st Qu.: 13.57   1st Qu.:12.50  
##  Median : 12081   Median :  55.0   Median : 35.00   Median :19.30  
##  Mean   : 17344   Mean   : 150.3   Mean   : 47.35   Mean   :20.88  
##  3rd Qu.: 23112   3rd Qu.: 190.0   3rd Qu.: 71.78   3rd Qu.:27.55  
##  Max.   :123124   Max.   :1100.0   Max.   :204.80   Max.   :57.50

The data is not normalised, so it is necesserily not really comparable.

5.2.2 Performing principal component analysis (PCA) on the raw (non-standardized) human data

We perform the PCA analysis for non-standardised data to see why it is important:

human_data_PCA_raw <- prcomp(human_data_)
summary(human_data_PCA_raw)

## Importance of components:
##                              PC1      PC2   PC3   PC4   PC5   PC6    PC7    PC8
## Standard deviation     1.821e+04 183.5476 24.84 11.23 3.696 1.539 0.1913 0.1568
## Proportion of Variance 9.999e-01   0.0001  0.00  0.00 0.000 0.000 0.0000 0.0000
## Cumulative Proportion  9.999e-01   1.0000  1.00  1.00 1.000 1.000 1.0000 1.0000

Most of the variance is included into the first component.

biplot(human_data_PCA_raw, choices = 1:2)

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length =
## arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length =
## arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length =
## arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length =
## arrow.len): zero-length arrow is of indeterminate angle and so skipped

The GNI_C is mostly contributing to the PC1. That is because GNI_C has values up to 123124, which is a lot more than any other variable. The nations with high GNI_C are to the left.

5.2.3 Standardising the variables in the human data and repeat the above analysis

As we saw earlier the data is to dissimilar to use for PCA, hence scaling:

human_data_scaled <- scale(human_data_)
summary(human_data_scaled)

##   FM_2ED_Ratio     FM_LFPR_Ratio         EduExp          LifeExp_B      
##  Min.   :-2.8446   Min.   :-2.5987   Min.   :-2.7620   Min.   :-2.7457  
##  1st Qu.:-0.5220   1st Qu.:-0.5286   1st Qu.:-0.6816   1st Qu.:-0.6272  
##  Median : 0.3476   Median : 0.2408   Median : 0.1131   Median : 0.2937  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.6123   3rd Qu.: 0.7484   3rd Qu.: 0.7381   3rd Qu.: 0.6524  
##  Max.   : 2.7134   Max.   : 1.6796   Max.   : 2.5239   Max.   : 1.4487  
##      GNI_C            MMRatio           AdBRate             RP_p        
##  Min.   :-0.9206   Min.   :-0.7127   Min.   :-1.1509   Min.   :-1.8531  
##  1st Qu.:-0.7057   1st Qu.:-0.6554   1st Qu.:-0.8315   1st Qu.:-0.7435  
##  Median :-0.2891   Median :-0.4549   Median :-0.3041   Median :-0.1398  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.3167   3rd Qu.: 0.1896   3rd Qu.: 0.6012   3rd Qu.: 0.5925  
##  Max.   : 5.8094   Max.   : 4.5338   Max.   : 3.8760   Max.   : 3.2512

human_data_PCA_scaled <- prcomp(human_data_scaled)
summary(human_data_PCA_scaled)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.0705 1.1430 0.87523 0.77704 0.66349 0.53451 0.45617
## Proportion of Variance 0.5359 0.1633 0.09575 0.07547 0.05503 0.03571 0.02601
## Cumulative Proportion  0.5359 0.6992 0.79496 0.87043 0.92546 0.96117 0.98718
##                            PC8
## Standard deviation     0.32021
## Proportion of Variance 0.01282
## Cumulative Proportion  1.00000

human_data_PCA_scaled

## Standard deviations (1, .., p=8):
## [1] 2.0705291 1.1430478 0.8752283 0.7770364 0.6634878 0.5345100 0.4561665
## [8] 0.3202124
## 
## Rotation (n x k) = (8 x 8):
##                       PC1         PC2          PC3         PC4        PC5
## FM_2ED_Ratio  -0.35624746  0.05544659 -0.258173727  0.61994816 -0.5923317
## FM_LFPR_Ratio  0.05140236  0.72511638 -0.574905230  0.04778188  0.2842048
## EduExp        -0.42844344  0.13794051 -0.068952381 -0.06994677  0.1577789
## LifeExp_B     -0.44425251 -0.02871318  0.112160198 -0.05130932  0.1525719
## GNI_C         -0.35074145  0.05123343 -0.198230209 -0.73273506 -0.4859385
## MMRatio        0.43708962  0.14765550 -0.123634905 -0.25658826 -0.1736624
## AdBRate        0.41035458  0.08931249  0.008796231  0.05375321 -0.4810245
## RP_p          -0.08404578  0.64720641  0.728586153  0.01511136 -0.1500635
##                       PC6          PC7         PC8
## FM_2ED_Ratio   0.19541082  0.057381277  0.16336822
## FM_LFPR_Ratio -0.03316792 -0.226716792 -0.07410761
## EduExp        -0.39597156  0.776784601 -0.05176554
## LifeExp_B     -0.42259100 -0.439295194  0.62590822
## GNI_C          0.11842241 -0.138379206 -0.16985528
## MMRatio        0.17148338  0.351654675  0.72304994
## AdBRate       -0.74968006 -0.078049728 -0.14549600
## RP_p           0.14102220  0.005543132 -0.02360083

The PC1 is still the most important component, but now it only explains around 54% of variance, PC2 captures around 16%, then PC3 less than 10% and all the other components display diminishing returns.

names_for_the_biplot <- round(summary(human_data_PCA_scaled)$importance[2, ] * 100, digits = 2)
names_for_the_biplot <- paste0(names(names_for_the_biplot), " (", names_for_the_biplot, "%)")

biplot(human_data_PCA_scaled, choices = 1:2, xlab = names_for_the_biplot[1], ylab = names_for_the_biplot[2], cex = c(0.5, 0.7))

These results make al ot more sense and also are in concordance with what I originally concluded form the first graphical analysis: 1. PC1 is negatively associated with GNI, Life expectancy, Education and the ratio of educated females in the - all correlted positivley, and increase in MMortality and adolecent births leads to increase in PC1. These values were correlated to each other. 2. PC2 increases with the number of female representatives and female percentage in labour force, which are correlated to each other. As mentioned above the differences are due to the GNI being a much larger influence without scaling. With scaling we can actually see what’s happening.

5.2.4 My interpretation of the data.

The PC1 is related to what is usually descirbed as development level of a nation. Developed nations have high income, better education and life expectancy, which developing conuthries and failed states have high adolescent pregnancies and mortality. It is widely believed, that female education is the main driving factor in the population reproducibility index decrease. The second PC shows that there is a connection between female representation in government and overall inclusion in work force, but it is not strictly connected ot the “PC1 - development level” meaning that in poorer countries females can also participate in the govrnment and workforce. Also in petrocraties of the Middle East, despite high GNI and associated parameters females have lower representation.

5.2.5 Exploring the tea data and analysis

library(FactoMineR)
tea_data <- read.csv("https://raw.githubusercontent.com/KimmoVehkalahti/Helsinki-Open-Data-Science/master/datasets/tea.csv", stringsAsFactors = TRUE)

# View(tea_data)

str(tea_data)

## 'data.frame':    300 obs. of  36 variables:
##  $ breakfast       : Factor w/ 2 levels "breakfast","Not.breakfast": 1 1 2 2 1 2 1 2 1 1 ...
##  $ tea.time        : Factor w/ 2 levels "Not.tea time",..: 1 1 2 1 1 1 2 2 2 1 ...
##  $ evening         : Factor w/ 2 levels "evening","Not.evening": 2 2 1 2 1 2 2 1 2 1 ...
##  $ lunch           : Factor w/ 2 levels "lunch","Not.lunch": 2 2 2 2 2 2 2 2 2 2 ...
##  $ dinner          : Factor w/ 2 levels "dinner","Not.dinner": 2 2 1 1 2 1 2 2 2 2 ...
##  $ always          : Factor w/ 2 levels "always","Not.always": 2 2 2 2 1 2 2 2 2 2 ...
##  $ home            : Factor w/ 2 levels "home","Not.home": 1 1 1 1 1 1 1 1 1 1 ...
##  $ work            : Factor w/ 2 levels "Not.work","work": 1 1 2 1 1 1 1 1 1 1 ...
##  $ tearoom         : Factor w/ 2 levels "Not.tearoom",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ friends         : Factor w/ 2 levels "friends","Not.friends": 2 2 1 2 2 2 1 2 2 2 ...
##  $ resto           : Factor w/ 2 levels "Not.resto","resto": 1 1 2 1 1 1 1 1 1 1 ...
##  $ pub             : Factor w/ 2 levels "Not.pub","pub": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Tea             : Factor w/ 3 levels "black","Earl Grey",..: 1 1 2 2 2 2 2 1 2 1 ...
##  $ How             : Factor w/ 4 levels "alone","lemon",..: 1 3 1 1 1 1 1 3 3 1 ...
##  $ sugar           : Factor w/ 2 levels "No.sugar","sugar": 2 1 1 2 1 1 1 1 1 1 ...
##  $ how             : Factor w/ 3 levels "tea bag","tea bag+unpackaged",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ where           : Factor w/ 3 levels "chain store",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ price           : Factor w/ 6 levels "p_branded","p_cheap",..: 4 6 6 6 6 3 6 6 5 5 ...
##  $ age             : int  39 45 47 23 48 21 37 36 40 37 ...
##  $ sex             : Factor w/ 2 levels "F","M": 2 1 1 2 2 2 2 1 2 2 ...
##  $ SPC             : Factor w/ 7 levels "employee","middle",..: 2 2 4 6 1 6 5 2 5 5 ...
##  $ Sport           : Factor w/ 2 levels "Not.sportsman",..: 2 2 2 1 2 2 2 2 2 1 ...
##  $ age_Q           : Factor w/ 5 levels "+60","15-24",..: 4 5 5 2 5 2 4 4 4 4 ...
##  $ frequency       : Factor w/ 4 levels "+2/day","1 to 2/week",..: 3 3 1 3 1 3 4 2 1 1 ...
##  $ escape.exoticism: Factor w/ 2 levels "escape-exoticism",..: 2 1 2 1 1 2 2 2 2 2 ...
##  $ spirituality    : Factor w/ 2 levels "Not.spirituality",..: 1 1 1 2 2 1 1 1 1 1 ...
##  $ healthy         : Factor w/ 2 levels "healthy","Not.healthy": 1 1 1 1 2 1 1 1 2 1 ...
##  $ diuretic        : Factor w/ 2 levels "diuretic","Not.diuretic": 2 1 1 2 1 2 2 2 2 1 ...
##  $ friendliness    : Factor w/ 2 levels "friendliness",..: 2 2 1 2 1 2 2 1 2 1 ...
##  $ iron.absorption : Factor w/ 2 levels "iron absorption",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ feminine        : Factor w/ 2 levels "feminine","Not.feminine": 2 2 2 2 2 2 2 1 2 2 ...
##  $ sophisticated   : Factor w/ 2 levels "Not.sophisticated",..: 1 1 1 2 1 1 1 2 2 1 ...
##  $ slimming        : Factor w/ 2 levels "No.slimming",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ exciting        : Factor w/ 2 levels "exciting","No.exciting": 2 1 2 2 2 2 2 2 2 2 ...
##  $ relaxing        : Factor w/ 2 levels "No.relaxing",..: 1 1 2 2 2 2 2 2 2 2 ...
##  $ effect.on.health: Factor w/ 2 levels "effect on health",..: 2 2 2 2 2 2 2 2 2 2 ...

Everything is a factor, apart from age. We can shoose the variables to look at. Seems more interesting to choose factual things, rather than emotional responses. The first 18 are the ones to choose from.

summary(tea_data)

##          breakfast           tea.time          evening          lunch    
##  breakfast    :144   Not.tea time:131   evening    :103   lunch    : 44  
##  Not.breakfast:156   tea time    :169   Not.evening:197   Not.lunch:256  
##                                                                          
##                                                                          
##                                                                          
##                                                                          
##                                                                          
##         dinner           always          home           work    
##  dinner    : 21   always    :103   home    :291   Not.work:213  
##  Not.dinner:279   Not.always:197   Not.home:  9   work    : 87  
##                                                                 
##                                                                 
##                                                                 
##                                                                 
##                                                                 
##         tearoom           friends          resto          pub     
##  Not.tearoom:242   friends    :196   Not.resto:221   Not.pub:237  
##  tearoom    : 58   Not.friends:104   resto    : 79   pub    : 63  
##                                                                   
##                                                                   
##                                                                   
##                                                                   
##                                                                   
##         Tea         How           sugar                     how     
##  black    : 74   alone:195   No.sugar:155   tea bag           :170  
##  Earl Grey:193   lemon: 33   sugar   :145   tea bag+unpackaged: 94  
##  green    : 33   milk : 63                  unpackaged        : 36  
##                  other:  9                                          
##                                                                     
##                                                                     
##                                                                     
##                   where                 price          age        sex    
##  chain store         :192   p_branded      : 95   Min.   :15.00   F:178  
##  chain store+tea shop: 78   p_cheap        :  7   1st Qu.:23.00   M:122  
##  tea shop            : 30   p_private label: 21   Median :32.00          
##                             p_unknown      : 12   Mean   :37.05          
##                             p_upscale      : 53   3rd Qu.:48.00          
##                             p_variable     :112   Max.   :90.00          
##                                                                          
##            SPC               Sport       age_Q          frequency  
##  employee    :59   Not.sportsman:121   +60  :38   +2/day     :127  
##  middle      :40   sportsman    :179   15-24:92   1 to 2/week: 44  
##  non-worker  :64                       25-34:69   1/day      : 95  
##  other worker:20                       35-44:40   3 to 6/week: 34  
##  senior      :35                       45-59:61                    
##  student     :70                                                   
##  workman     :12                                                   
##              escape.exoticism           spirituality        healthy   
##  escape-exoticism    :142     Not.spirituality:206   healthy    :210  
##  Not.escape-exoticism:158     spirituality    : 94   Not.healthy: 90  
##                                                                       
##                                                                       
##                                                                       
##                                                                       
##                                                                       
##          diuretic             friendliness            iron.absorption
##  diuretic    :174   friendliness    :242   iron absorption    : 31   
##  Not.diuretic:126   Not.friendliness: 58   Not.iron absorption:269   
##                                                                      
##                                                                      
##                                                                      
##                                                                      
##                                                                      
##          feminine             sophisticated        slimming          exciting  
##  feminine    :129   Not.sophisticated: 85   No.slimming:255   exciting   :116  
##  Not.feminine:171   sophisticated    :215   slimming   : 45   No.exciting:184  
##                                                                                
##                                                                                
##                                                                                
##                                                                                
##                                                                                
##         relaxing              effect.on.health
##  No.relaxing:113   effect on health   : 66    
##  relaxing   :187   No.effect on health:234    
##                                               
##                                               
##                                               
##                                               
## 

dim(tea_data)

## [1] 300  36

keep_tea_columns <- c("Tea", "How", "how", "sex",  "frequency","age_Q")



# select the 'keep_columns' to create a new dataset
tea_time <- dplyr::select(tea_data,  all_of(keep_tea_columns))

I choose theese variables to look at.

ggpairs(tea_time,progress = F) 

summary(tea_time)

##         Tea         How                      how      sex           frequency  
##  black    : 74   alone:195   tea bag           :170   F:178   +2/day     :127  
##  Earl Grey:193   lemon: 33   tea bag+unpackaged: 94   M:122   1 to 2/week: 44  
##  green    : 33   milk : 63   unpackaged        : 36           1/day      : 95  
##                  other:  9                                    3 to 6/week: 34  
##                                                                                
##    age_Q   
##  +60  :38  
##  15-24:92  
##  25-34:69  
##  35-44:40  
##  45-59:61

People like Earl Grey. There are prefered methods of drinking tea.

mca <- MCA(tea_time)

summary(mca)

## 
## Call:
## MCA(X = tea_time) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5   Dim.6   Dim.7
## Variance               0.262   0.241   0.210   0.196   0.194   0.184   0.174
## % of var.             10.460   9.638   8.416   7.858   7.744   7.375   6.950
## Cumulative % of var.  10.460  20.098  28.514  36.372  44.115  51.490  58.440
##                        Dim.8   Dim.9  Dim.10  Dim.11  Dim.12  Dim.13  Dim.14
## Variance               0.169   0.152   0.143   0.142   0.122   0.112   0.108
## % of var.              6.777   6.069   5.723   5.674   4.894   4.493   4.324
## Cumulative % of var.  65.217  71.287  77.010  82.683  87.578  92.070  96.395
##                       Dim.15
## Variance               0.090
## % of var.              3.605
## Cumulative % of var. 100.000
## 
## Individuals (the 10 first)
##                       Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
## 1                  |  0.522  0.347  0.113 |  0.206  0.059  0.018 | -0.471
## 2                  |  0.406  0.210  0.069 | -0.092  0.012  0.004 | -0.982
## 3                  | -0.227  0.066  0.040 | -0.383  0.203  0.112 | -0.064
## 4                  | -0.534  0.363  0.221 |  0.338  0.158  0.088 | -0.377
## 5                  | -0.029  0.001  0.001 |  0.058  0.005  0.002 | -0.140
## 6                  | -0.534  0.363  0.221 |  0.338  0.158  0.088 | -0.377
## 7                  |  0.113  0.016  0.004 |  0.120  0.020  0.005 |  0.134
## 8                  |  0.486  0.302  0.069 | -0.258  0.092  0.019 | -0.746
## 9                  |  0.165  0.035  0.010 | -0.072  0.007  0.002 | -0.257
## 10                 |  0.576  0.422  0.132 | -0.267  0.099  0.028 |  0.273
##                       ctr   cos2  
## 1                   0.352  0.092 |
## 2                   1.527  0.403 |
## 3                   0.007  0.003 |
## 4                   0.225  0.110 |
## 5                   0.031  0.014 |
## 6                   0.225  0.110 |
## 7                   0.028  0.006 |
## 8                   0.881  0.162 |
## 9                   0.105  0.025 |
## 10                  0.118  0.030 |
## 
## Categories (the 10 first)
##                        Dim.1     ctr    cos2  v.test     Dim.2     ctr    cos2
## black              |   1.178  21.809   0.454  11.654 |  -0.433   3.192   0.061
## Earl Grey          |  -0.530  11.509   0.506 -12.304 |   0.021   0.020   0.001
## green              |   0.457   1.467   0.026   2.780 |   0.845   5.427   0.088
## alone              |  -0.198   1.620   0.073  -4.659 |   0.014   0.009   0.000
## lemon              |   0.169   0.200   0.004   1.027 |   0.272   0.563   0.009
## milk               |   0.251   0.841   0.017   2.235 |   0.136   0.270   0.005
## other              |   1.910   6.974   0.113   5.808 |  -2.256  10.559   0.157
## tea bag            |  -0.146   0.771   0.028  -2.889 |   0.007   0.002   0.000
## tea bag+unpackaged |  -0.194   0.749   0.017  -2.261 |  -0.198   0.849   0.018
## unpackaged         |   1.195  10.929   0.195   7.633 |   0.485   1.953   0.032
##                     v.test     Dim.3     ctr    cos2  v.test  
## black               -4.280 |  -0.072   0.102   0.002  -0.715 |
## Earl Grey            0.498 |  -0.004   0.001   0.000  -0.094 |
## green                5.134 |   0.186   0.301   0.004   1.129 |
## alone                0.330 |   0.328   5.524   0.199   7.718 |
## lemon                1.654 |   0.582   2.951   0.042   3.538 |
## milk                 1.216 |  -1.200  23.938   0.383 -10.694 |
## other               -6.859 |  -0.834   1.651   0.021  -2.535 |
## tea bag              0.133 |  -0.563  14.221   0.414 -11.130 |
## tea bag+unpackaged  -2.312 |   0.555   7.642   0.140   6.481 |
## unpackaged           3.097 |   1.209  13.897   0.199   7.720 |
## 
## Categorical variables (eta2)
##                      Dim.1 Dim.2 Dim.3  
## Tea                | 0.546 0.125 0.005 |
## How                | 0.151 0.165 0.430 |
## how                | 0.195 0.041 0.451 |
## sex                | 0.089 0.407 0.011 |
## frequency          | 0.021 0.286 0.270 |
## age_Q              | 0.567 0.422 0.096 |

plot(mca, graph.type="ggplot", invisible=c("var","quali.sup","quanti.sup"),cex=0.8)

plot(mca, graph.type="ggplot", invisible=c("ind"), cex=0.8, habillage="quali")

Looking at the plot we can assume that there is a correlation between the age and which tea the person prefers:

• 15-24 prefer earl gray
• oldest prefer black tea.
• 25-35 prefer green tea
• Drinking tea straight is a younger quality, older prefer tea with something.
• Middle age is associated with green tea bags daily.